Contaminantes y zonas verdes
Fernando Martínez Gómez
library(dplyr, warn.conflicts = FALSE)
library(tidyr, warn.conflicts = FALSE)
library(arrow, warn.conflicts = FALSE)
library(leaflet, warn.conflicts = FALSE)
library(sf, warn.conflicts = FALSE)## Linking to GEOS 3.12.2, GDAL 3.8.4, PROJ 9.4.0; sf_use_s2() is TRUE## WARNING: different compile-time and runtime versions for GEOS found:## Linked against: 3.12.2-CAPI-1.18.2 compiled against: 3.12.1-CAPI-1.18.1## It is probably a good idea to reinstall sf (and maybe lwgeom too)library(RColorBrewer, warn.conflicts = FALSE)
library(htmltools, warn.conflicts = FALSE)
library(ggplot2, warn.conflicts = FALSE)
library(scales, warn.conflicts = FALSE)
library(lubridate)##
## Adjuntando el paquete: 'lubridate'## The following object is masked from 'package:arrow':
##
## duration## The following objects are masked from 'package:base':
##
## date, intersect, setdiff, unionlibrary(knitr, warn.conflicts = FALSE)z_verdes = read_parquet("../../../data/cleaned/zonas_verdes/espacios_verdes_limpio.parquet")1 Cruze de Zonas Verdes y Contaminantes
1.1 Pivotación de Zonas Verdes
Vamos a agrupar por distritos municipales y sumar las áreas de las zonas verdes.
Antes de hacer eso, hay un distrito ficticio que es el cauce del río Turia. Como tenemos la localización de las zonas, podemos reclasificarlo en el distrito real que corresponda, dependiendo de las coordenadas.
Como solo hay 15 zonas verdes en el cauce del Turia, podemos hacer una reclasificación manual.
z_verdes_2 <- z_verdes %>%
mutate(
Distrito_Municipal = case_when(
OBJECTID == 1257 ~ "PLA DEL REIAL",
OBJECTID == 1570 ~ "CAMINS AL GRAU",
OBJECTID == 1728 ~ "CAMINS AL GRAU",
OBJECTID == 1214 ~ "SAIDIA",
OBJECTID == 1233 ~ "OLIVERETA",
OBJECTID == 1305 ~ "CIUTAT VELLA",
OBJECTID == 1494 ~ "POBLATS MARITIMS",
OBJECTID == 1493 ~ "PLA DEL REIAL",
OBJECTID == 1709 ~ "QUATRE CARRERES",
OBJECTID == 1195 ~ "EXTRAMURS",
OBJECTID == 1288 ~ "EXTRAMURS",
OBJECTID == 1295 ~ "OLIVERETA",
OBJECTID == 1471 ~ "SAIDIA",
OBJECTID == 1508 ~ "SAIDIA",
OBJECTID == 1799 ~ "PLA DEL REIAL",
TRUE ~ Distrito_Municipal # Mantener el valor original si no coincide con ningún OBJECTID
)
)z_verdes_distritos <- z_verdes_2 %>%
group_by(Distrito_Municipal) %>%
summarise(area_total = sum(Area)) %>%
ungroup()
z_verdes_distritosdistritos_sf <- st_read("../../../data/cleaned/geoespacial/districtes-distritos.geojson", quiet = TRUE)
distritos_sfHay que cambiar el nombre de algunos distritos municipales para que coincida con el nombre en distritos_sf.
z_verdes_distritos_2 <- z_verdes_distritos %>%
mutate(
Distrito_Municipal = dplyr::recode(
Distrito_Municipal,
"POBLATS DE L`OEST" = "POBLATS DE L'OEST",
"OLIVERETA" = "L'OLIVERETA",
"SAIDIA" = "LA SAIDIA",
"PLA DEL REIAL" = "EL PLA DEL REAL",
"L´EIXAMPLA" = "L'EIXAMPLE",
"POBLES DEL SUD" = "POBLATS DEL SUD"
)
)También hay que dejar una sola fila de “POBLATS DEL NORD” en distritos_sf, porque hay tres filas con el mismo nombre.
distritos_sf <- distritos_sf %>%
filter(!(objectid %in% c(217, 199)))distritos_sf <- distritos_sf %>%
left_join(z_verdes_distritos_2, by = c("nombre" = "Distrito_Municipal"))2 Mapa zonas verdes
pal <- colorNumeric(
palette = "YlGnBu",
domain = distritos_sf$area_total,
na.color = "transparent"
)leaflet(distritos_sf) %>%
addProviderTiles(providers$CartoDB.Positron) %>%
addPolygons(
fillColor = ~pal(area_total),
weight = 1,
opacity = 1,
color = "white",
dashArray = "3",
fillOpacity = 0.7,
highlightOptions = highlightOptions(
weight = 3,
color = "#666",
dashArray = "",
fillOpacity = 0.9,
bringToFront = TRUE
),
label = ~paste0(
"<strong>", nombre, "</strong><br/>",
"Superficie parques: ",
formatC(area_total, big.mark = ".", format = "d"), " m²"
) %>% lapply(htmltools::HTML),
labelOptions = labelOptions(
style = list("font-weight" = "normal", padding = "3px 8px"),
textsize = "13px",
direction= "auto"
)
) %>%
addLegend(
pal = pal,
values = ~area_total,
opacity = 0.7,
title = "m² parques",
position = "bottomright"
)3 Mapa contaminantes
contaminantes = read_parquet("../../../data/cleaned/calidad_aire/conjunto_estaciones_limpio.parquet")
contaminantescontaminantes = contaminantes %>%
mutate(
FechaHora = as.POSIXct(
FechaHora,
format = "%Y-%m-%d %H:%M:%S",
tz = "Europe/Madrid"
)
)
codigos = read.csv("../../../data/cleaned/geoespacial/codigos.csv")contaminantes_23 <- contaminantes %>%
filter(year(FechaHora) == 2023) contaminantes_23 <- contaminantes_23 %>%
mutate(
Origen = dplyr::recode(
Origen,
"Massanassa_UM" = "Massanassa",
)
)zonas_sel <- c("Massanassa", "Quart de Poblet", "Sedavi UM",
"Torrent-El Vedat", "València - Av. França",
"València - Molí del Sol", "València - Pista de Silla",
"València - Politècnic", "València Port llit antic Túria",
"València Port Moll Trans. Ponent")contaminantes_23_cp = contaminantes_23 %>%
filter(Origen %in% zonas_sel) %>%
left_join(codigos, by = c("Origen" = "Origen")) %>%
group_by(COD_POSTAL) %>%
summarise(
PM10 = mean(PM10, na.rm = TRUE),
PM2.5 = mean(PM2.5, na.rm = TRUE),
SO2 = mean(SO2, na.rm = TRUE),
NO = mean(NO, na.rm = TRUE),
NO2 = mean(NO2, na.rm = TRUE),
NOx = mean(NOx, na.rm = TRUE),
O3 = mean(O3, na.rm = TRUE),
.groups = "drop"
)# Eliminar filas cuyo COD_POSTAL sea 46901 y 46930
contaminantes_23_dist = contaminantes_23_cp %>%
filter(!(COD_POSTAL %in% c(46901, 46930))) %>%
mutate(
distrito = case_when(
COD_POSTAL == 46001 ~ "CIUTAT VELLA",
COD_POSTAL %in% c(46015, 46023) ~ "CAMPANAR",
COD_POSTAL == 46022 ~ "AYORA",
COD_POSTAL == 46024 ~ "QUATRE CARRERES"
)
)Ya podemos crear el choroplete de contaminantes por distrito.
contaminantes_23_dist = distritos_sf %>%
left_join(contaminantes_23_dist, by = c("nombre" = "distrito"))Mapa de NO2
Mapa de NOx
Mapa de NO
Mapa de SO2
Mapa de PM10
Mapa de PM2.5
Mapa de O3
zverdes_contaminantes <- distritos_sf %>%
left_join(
contaminantes_23_dist %>% st_set_geometry(NULL), # quitamos geometría duplicada
by = "nombre"
)## Warning in sf_column %in% names(g): Detected an unexpected many-to-many relationship between `x` and `y`.
## ℹ Row 6 of `x` matches multiple rows in `y`.
## ℹ Row 6 of `y` matches multiple rows in `x`.
## ℹ If a many-to-many relationship is expected, set `relationship =
## "many-to-many"` to silence this warning.4 Correlación entre áreas verdes y contaminantes
df_vals <- zverdes_contaminantes %>%
st_set_geometry(NULL) %>%
select(area_total.x, NO2, NOx, NO, SO2, PM10, `PM2.5`, O3)
# 4. Matriz de correlaciones
cor_mat <- cor(df_vals, use = "pairwise.complete.obs")
# 5. Formato “largo” y filtrado para area_total.x
cor_df <- as.data.frame(cor_mat) %>%
tibble::rownames_to_column("var1") %>%
pivot_longer(
cols = -var1,
names_to = "var2",
values_to = "corr"
) %>%
# Solo mantener las correlaciones de area_total.x
filter(var1 == "area_total.x" & var2 != "area_total.x")
# 6. Gráfico de barras de correlaciones
ggplot(cor_df, aes(x = var2, y = corr, fill = corr)) +
geom_col() +
scale_fill_gradient2(
low = "blue",
mid = "white",
high = "red",
midpoint = 0,
limits = c(-1, 1),
name = "ρ"
) +
geom_text(aes(label = sprintf("%.2f", corr)), vjust = ifelse(cor_df$corr>0, -0.5, 1.5)) +
theme_minimal(base_size = 14) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
axis.title = element_blank()
) +
labs(title = "Correlación de area_total.x con contaminantes")
df_vals <- zverdes_contaminantes %>%
st_set_geometry(NULL) %>%
select(area_total.x, NO2, NOx, NO, SO2, PM10, `PM2.5`, O3)
# 2. Definimos el vector de contaminantes
contaminantes <- c("NO2","NOx","NO","SO2","PM10","PM2.5","O3")
# 3. Correlaciones y tests en bucle
resumen_cor <- lapply(contaminantes, function(var) {
test <- cor.test(df_vals$area_total.x, df_vals[[var]], method = "pearson")
# tidy() de broom convierte el resultado en data.frame
broom::tidy(test) %>%
select(estimate, p.value, parameter) %>%
mutate(
contaminante = var,
n = parameter + 2,
signif = if_else(p.value < 0.05, "p<0.05", "n.s.")
) %>%
rename(
r = estimate,
df = parameter
)
}) %>%
bind_rows() %>%
select(contaminante, n, df, r, p.value, signif) %>%
arrange(desc(abs(r)))
print(resumen_cor)## # A tibble: 7 × 6
## contaminante n df r p.value signif
## <chr> <dbl> <int> <dbl> <dbl> <chr>
## 1 PM2.5 4 2 0.837 0.163 n.s.
## 2 NO2 4 2 0.629 0.371 n.s.
## 3 PM10 4 2 0.474 0.526 n.s.
## 4 SO2 4 2 -0.387 0.613 n.s.
## 5 NOx 4 2 0.252 0.748 n.s.
## 6 NO 4 2 0.170 0.830 n.s.
## 7 O3 4 2 -0.0244 0.976 n.s.